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Many small moonlets, creating propeller structures, have been found in Saturn's rings by the Cassini spacecraft. We study the dynam- 
ical evolution of such 20-50m sized bodies which are embedded in Saturn's rings. We estimate the importance of various interaction 
processes with the ring particles on the moonlet's eccentricity and semi-major axis analytically. For low ring surface densities, the 
main effects on the evolution of the eccentricity and the semi-major axis are found to be due to collisions and the gravitational interac- 
tion with particles in the vicinity of the moonlet. For large surface densities, the gravitational interaction with self-gravitating wakes 
becomes important. 

We also perform realistic three dimensional, collisional N-body simulations with up to a quarter of a million particles. A new set 
of pseudo shear periodic boundary conditions is used which reduces the computational costs by an order of magnitude compared to 
previous studies. Our analytic estimates are confirmed to within a factor of two. 

On short timescales the evolution is always dominated by stochastic effects caused by collisions and gravitational interaction with self- 
gravitating ring particles. These result in a random walk of the moonlet's semi-major axis. The eccentricity of the moonlet quickly 
reaches an equilibrium value due to collisional damping. The average change in semi-major axis of the moonlet after 100 orbital 
periods is 10- 100m. This translates to an offset in the azimuthal direction of several hundred kilometres. We expect that such a shift 
is easily observable. 
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1. Introduction 

Small bodies, so called moonlets, which are embedded in 
Saturn's rings, can create propeller shaped structures due to 
their disturbance of the rings. These have been predicted 
both analytically and n umerically llSpahn and Sremcevid |2000L 
ISremcevic etail l2002l ISeifi et all T2005H . Only recently, they 



have b een observed by the Cassini s pacecraft in both the A and 
B ring llTiscareno et al.Ll2006ll2008ll . 

These 20 m - 100 m sized bodies can migrate within the 
rings, similar to proto-planets which migrate in a proto-stellar 
disc. Depending on the disc properties and the moonlet size, this 
can happen in either a smooth or in a stochastic (random walk) 
fashion. We refer to those migration regimes as type I and type 
IV, respec tively, in analogy t o the terminology in disc-planet in- 
teractions. jCrid^etal] Il2010ll showed that there is a laminar type 
I regime that might be important on very long timescales. This 
migration is qualitatively different to the migration in a pressure 
supported gas disc. However, the migration of moonlets in the 
A ring is generally dominated by type IV migration, at least on 
short timescales. 

In this paper, we study the type IV migration regime in de- 
tail. The full mutual interaction of the ring particles with the 
moonlet and its consequent induced motion are considered both 
analytically and numerically. 

We first review the basic equations governing the moonlet 
and the ring particles in a shearing box approximation in Sect. [2] 
Then, in Sect.[3j we estimate the eccentricity damping timescale 
due to ring particles colliding with the moonlet and ring particles 
on horseshoe orbits as well as the effect of particles on circulat- 
ing orbits. In Sect. [4] we estimate the excitation of the moonlet 



eccentricity caused by stochastic particle collisions and gravita- 
tional interactions with ring particles. This enables us to derive 
an analytic estimate of the equilibrium eccentricity. 

In Sect. |5] we discuss and evaluate processes, such as colli- 
sions and gravitational interactions with ring particles and self- 
gravitating clumps, that lead to a random walk in the semi-major 
axis of the moonlet. 

We describe our numerical code and the initial conditions 
used in Sect. [6] We perform realistic three dimensional N-body 
simulations of the ring system and the moonlet, taking into ac- 
count a moonlet with finite size, a size distribution of ring par- 
ticles, self-gravity and collisions. The results are presented in 
Sect. [7] All analytic estimates are confirmed both in terms of 
qualitative trends and quantitatively to within a factor of about 
two in all simulations that we performed. We also discuss the 
long term evolution of the longitude of the moonlet and its ob- 
servability, before summarising our results in Sec. [8] 



2. Basic equations governing the moonlet and ring 
particles 

We adopt a local right handed Cartesian coordinate system with 
its origin being in circular Keplerian orbit with semi-major axis 
a and rotating uniformly with angular velocity Q. This orbit co- 
incides with that of the moonlet when it is assumed to be un- 
perturbed by ring particles. The x axis coincides with the line 
joining the central object of mass M p and the origin. The unit 
vector in the x direction, e c , points away from the central ob- 
ject. The unit vector in the y direction e v points in the direction 
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Fig. 1. Trajectories of ring particles in a frame centred on the moonlet. Particles accumulate near the moonlet and fill it's Roche 
lobe. Particles on trajectories labelled a) are on circulating orbits. Particles on trajectories labelled b) collide with other particles in 
the moonlet's vicinity. Particles on trajectories labelled (c) collide with the moonlet directly. Particles on trajectories labelled (d) are 
on horseshoe orbits. Particles on trajectories labelled (e) leave the vicinity of the moonlet. 



of rotation and the unit vector in the z direction, e ; points in the 
vertical direction being normal to the disc mid-plane. 

In general, we shall consider a ring particle of mass mi in- 
teracting with the moonlet which has a much larger mass m 2 . A 
sketch of three possible types of particle trajectory in the vicin- 
ity of the moonlet is shown in Fig.Q] These correspond to three 
distinct regimes, a) denoting circulating orbits, b) and c) denot- 
ing orbits that result in a collision with particles in the vicinity 
of the moonlet and directly with the moonlet respectively, and 
d) denoting horseshoe orbits. All these types of trajectory occur 
for particles that are initially in circular orbits, both interior and 
exterior to the moonlet, when at large distances from it. 

Approximating the gravitational force due to the central ob- 
ject by its first order Taylor expansion about the origin leads to 
Hill's equation, governing the motion of a particle of mass m\ of 
the form 



plus the circular Keplerian velocity corresponding to the orbital 
frequency Q. of the origin. Thus 



= |i, y + + a Q Z 



(4) 



All eccentricities considered here are very small (~ 10~ 8 - 1CT 7 ) 
so that the difference between the quantities defined through use 
of Eqs.|2]and[3]is negligible, allowing us to adopt £ as a measure 
of the eccentricity throughout this paper. 

Let us define another quantity, 3K, that is also conserved for 
non interacting particle motion, Y = 0, which is the x coordinate 
of the centre of the epicyclic motion and is given by 



M = 2Cr 1 y + 4x. 



(5) 



r = -2Qe z x r + 3Q 2 (r • eje.t - V*r/mi 



(1) 



We identify a change in J[ as a change in the semi-major axis a 
of the particle, again, under the assumption that the eccentricity 
is small. 



where r = (x,y,z) is position of a particle with mass ni\ and *F 
is the gravitational potential acting on the particle due to other 
objects of interest such as the moonlet. 

The square amplitude of the epicyclic motion £ 2 can be de- 
fined through 



£T 2 x 2 + (2QT l y + 3x) 2 . 



(2) 



Note that neither the eccentricity e, nor the semi-major axis a 
are formally defined in the local coordinate system. However, 
in the absence of interaction with other masses QP = 0), £ is 
conserved, and up to first order in the eccentricity we may make 
the identification £ = ea. We recall the classical definition of 
the eccentricity e in a coordinate system centred on the central 
object 



w x (s x w) 



Q 2 a 3 



(3) 



Here, s is the position vector (a, 0, 0) and w is the velocity vector 
of the particle relative to the mean shear associated with circu- 
lar Keplerian motion as viewed in the local coordinate system 



2. 1 . Two interacting particles 

We now consider the motion of two gravitationally interact- 
ing particles with position vectors ri = Qci.yi.Zi) and r 2 = 
{x2,y 2 ,z 2 ). Their corresponding masses are m\ and m%, respec- 
tively. The governing equations of motion are 



x x = -2Qe ; x tj + 3Q 2 (r, ■ e,)e. v - V r ,>F 12 /mi 
f 2 = -2Qe ; x r 2 + 3Q 2 (r 2 • e x )e x + V r y l2 /m 2 , 



(6) 
(7) 



where the interaction gravitational potential is ¥12 = 
-Gmxmil \y\ - Y2\- The position vector of the centre of mass of 
the two particles is given by 



m\t\ + m 2 r 2 



(8) 



mi + m 2 

The vector f also obeys Eq.Q]with *P = 0, which applies to an 
isolated particle. This is because the interaction potential does 
not affect the motion of the centre of mass. We also find it useful 
to define the vector 



= (tr'i,-, 2fr 1 j / + 3xd i = 1,2. 



(9) 
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Then, consistently with our earlier definition of &, we have <5, = 
|£,|. The amplitude of the epicyclic motion of the centre of mass 
S is given by 

(10) 



(nt\ + mi) 2 & 2 



m\s\ + m\&\ + 7m\m%£i\&i cos(0i2). 



Here (f>n is the angle between S\ and Si. It is important to note 
that & is conserved even if the two particles approach each other 
and become bound. This is as long as frictional forces are in- 
ternal to the two particle system and do not affect the centre of 
mass motion. 

3. Effects leading to damping of the eccentricity of 
the moonlet 

We begin by estimating the moonlet eccentricity damping rate 
associated with ring particles that either collide directly with the 
moonlet, particles in its vicinity, or only interact gravitationally 
with the moonlet. 



3. 1 . The contribution due to particles impacting the moonlet 

Particles impacting the moonlet in an eccentric orbit exchange 
momentum with it. Let us assume that a ring particle, identified 
with mi has zero epicyclic amplitude, so that S\ = far away 
from the moonlet. The moonlet is identified with mi and has 
an initial epicyclic amplitude Si. The epicyclic amplitude of the 
centre of mass is therefore 



rti\ + mi \ mi) 



(11) 



where we have assumed that ni\ <K mi. 

The moonlet is assumed to be in a steady state in which there 
is no net accretion of ring particles. Therefore, for every particle 
that either collides directly with the moonlet or nearby particles 
bound to it (and so itself becomes temporarily bound to it), one 
particle must also escape from the moonlet. This happens pri- 
marily through slow leakage from locations close to the Li and 
L2 points such that most particles escape from the moonlet with 
almost zero velocity (as viewed from the centre of mass frame) 
and so do not change its orbital eccentricity. However, after an 
impacting particle becomes bound to the moonlet, conservation 
of the epicyclic amplitude associated with the centre of mass 
motion together with Eq.Ql] imply that each impacting particle 
will reduce the eccentricity by a factor 1 - m\jmi. 

It is now an easy task to estimate the eccentricity damping 
timescale by determining the number of particle collisions per 
time unit with the moonlet or particles bound to it. To do that, a 
smooth window function Wb+ C (b) is used, being unity for impact 
parameters b that always result in an impact with the moonlet 
or particles nearby that are bound to it, being zero for impact 
parameters that never result in an impact. 

The number of particles impacting the moonlet per time unit, 
dN/dt, is obtained by integrating over the impact parameter with 
the result that 



— = — -i.n\b\w b+c (b)db. 

at mi J_ M 2 



(12) 



We note that allowing b to be negative enables impacts from both 
sides of the moonlet to be taken into account. 

Therefore, after using Eq.QT|we find that the rate of change 
of the moonlet's eccentricity e%, or equivalently of it's amplitude 
of epicyclic motion Si, is given by 

dSo Si So r°° 3 

= -— -2Q|ft W b+C (b)db, (13) 
m 2 J_ M 2 



dt 



where t<, collisions defines the circularisation time arising from col- 
lisions with the moonlet. We remark that the natural unit for b is 
the Hill radius of the moonlet, r# = (mi/(3M p )) l/3 a, so that the 
dimensional scaling for T e>C oiiisions is given by 



Collisions KGIf H lQ ' 



(14) 



which we find to also apply to all the processes for modifying 
the moonlet's eccentricity discussed below. If we assume that 
WW(|&|) can be approximated by a box function, being unity in 
the interval [1.5r#, 2.5 r#] and zero elsewhere, we get 



Collisions = 2.0 GZr^n- 1 . 



(15) 



3.2. Eccentricity damping due the interaction of the moonlet 
with particles on horseshoe orbits 

The eccentricity of the moonlet manifests itself in a small oscil- 
lation of the moonlet about the origin. Primarily ring particles 
on horseshoe orbits will respond to that oscillation and damp it. 
This is because only those particles on horseshoe orbits spend 
enough time in the vicinity of the moonlet, i.e. many epicyclic 
periods. 

In appendix [A] we calculate the amplitude of epicyclic mo- 
tion &\f (or equivalently the eccentricity e\f) that is induced in 
a single ring particle in a horseshoe orbit undergoing a close ap- 
proach to a moonlet which is assumed to be in an eccentric orbit. 
In order to calculate the circularisation time, we have to consider 
all relevant impact parameters. We begin by noting that each par- 
ticle encounter with the moonlet is conservative and is such that 
for each particle, the Jacobi constant, applicable when the moon- 
let is in circular orbit, is increased by an amount m\il 2 S\^2 by 
the action of the perturbing force, associated with the eccentric- 
ity of the moonlet, as the particle passes by. Because the Jacobi 
constant, or energy in the rotating frame, for the moonlet and the 
particle together is conserved, the square of the epicyclic ampli- 
tude associated with the moonlet alone changes by S\^m\jmi. 
Accordingly the change in the amplitude of epicyclic motion 
of the moonlet Si, consequent on inducing the amplitude of 
epicyclic motion S\f in the horseshoe particle, is given by 



A6 2 = _^i £ 2 

mi 1 



(16) 



Note that this is different compared to Eq.QT] Here, we are deal- 
ing with a second order effect. First, the eccentric moonlet ex- 
cites eccentricity in a ring particle. Second, because the total 
epicyclic motion is conserved, the epicyclic motion of the moon- 
let is reduced. 

Integrating over the impact parameters associated with ring 
particles and taking into account particles streaming by the 
moonlet from both directions by allowing negative impact pa- 
rameters gives 



dS\ 
dt 



horseshoe 



U —1 



°° 3Y£\p.\b\W d {b)db 
2mo 



2S\ 



' e, horse shoe 



-,(17) 



where T e _horseshoe is the circularisation time and, as above, we 
have inserted a window function, which is unity on impact pa- 
rameters that lead to horseshoe orbits, otherwise being zero. 
Using £1 / given by Eq. IA.20l we obtain 



~Te, collisions 



*~e,horseshoe 



128 



CLr 2 H Q.^ 

mi 



jy /3 W d ((2770,) 1/6 ) d V . (18) 

CO 
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For a sharp cutoff of Wd(b) at b m = 1.5rn we find the value of 
the integral in the above to be 2.84. Thus, in this case we get 



:Lseshoe = 0.13 GEr^Cr 1 . 



(19) 



However, note that this value is sensitive to the value of b m 
adopted. For b m = 1 25r H , t c is a factor of 4.25 larger. 

3.3. The effects of circulating particles 

The effect of the response of circulating particles to the gravi- 
tational perturbation of the moonlet on the moonlet's eccentric- 
ity can be estimated from the work of Goldreich and Tremaind 
1 198011 and lGoldreich and Tremainel lll982ll . These authors con- 
sidered a ring separated from a satellite such that co-orbital ef- 
fects were not considered. Thus, their expressions may be ap- 
plied to estimate effects due to circulating particles. However, 
we exclude their corotation torques as they are determined by 
the ring edges and are absent in a local model with uniform az- 
imuthally averaged surface density. Equivalently, one may sim- 
ply assume that the corotation torques are saturated. When this 
is done only Lindblad torques act on the moonlet. These tend to 
excite the moonlet's eccentricity rather than damp it. 

We replace the ring mass in Eq. 70 of 
iGoldreich and Tremaind Ill982ll by the integral over the 
impact parameter and switch to our notation to obtain 



d&j 
It 



= 9.55 m 2 EG 2 fr 



-5 c-2 



&i W a (b) db, 



(20) 



where W a (b) is the appropriate window function for circulating 
particles. Assuming a sharp cutoff of W a (|fe|) at b m , we can eval- 
uate the integral in Eq. [20]to get 



1 d&i 

~&L IT 



1 



*e,circ 



0.183 CLrT, 1 QT 1 , 



(21) 



where we have adopted b m = 2.5r#, consistently with the simu- 
lation results presented below. We see that, although T e ^ rc < 
corresponds to growth rather than damping of the eccentricity, it 
scales in the same manner as the circularisation times in Sects. 
Oand[32]scale (see Eq.fBl. 

This timescale and the timescale associated with particles 
on horseshoe orbits, T Cj horseshoe, are significantly larger than the 
timescale associated with collisions, given in Eq. [15] We can 
therefore ignore those effects for most of the discussion in this 
paper. 

4. Processes leading to the excitation of the 
eccentricity of the moonlet 

In Sec. [3] we assumed that the moonlet had a small eccentric- 
ity but neglected the initial eccentricity of the impacting ring 
particles. When this is included, collisions and gravitational in- 
teractions of ring particles with the moonlet may also excite its 
eccentricity. 

4.1. Collisional eccentricity excitation 

To see this, assume that the moonlet initially has zero eccentric- 
ity, or equivalently no epicyclic motion, but the ring particles do. 

As above we consider the conservation of the epicyclic mo- 
tion of the centre of mass in order to connect the amplitude of the 
final epicyclic motion of the combined moonlet and ring parti- 
cle to the initial amplitude of the ring particle's epicyclic motion. 



This gives the epicyclic amplitude of the centre of mass after one 
impact as 



£ = 



m.2 \ mo 



(22) 



Assuming that successive collisions are uncorrected and occur 
stochastically with the mean time between consecutive encoun- 
ters being r ce , the evolution of <5 is governed by the equation 

\2 



dt 



collisions 
-1 _ 



2\_-l 



(23) 



where r~ e = dN/dt can be expressed in terms of the surface 
density and an impact window function Wb+ C (b) (see Eq. [121 . 
The quantity (S 2 ) is the mean square value of Si for the ring 
particles. Using Eq.[2] this may be written in terms of the mean 
squares of the components of the velocity dispersion relative to 
the background shear, in the form 

(&]) = Q~ 2 <(i 2 +4(y l + 3Q Xl /2) 2 >. (24) 

We find in numerical simulations (ei) ~ 10~ 6 for almost all ring 
parameters. This value is mai nly determined by the coef ficient 
of restitution [see e.g. Fig 4 in lMorishima and Salol I2006T1 . 

4.2. Stochastic excitation due to circulating particles 

Ring particles that are on circulating orbits, such as that given 
by path (a) in Fig. Q] exchange energy and angular momen- 
tum with the moonlet and th erefore change its eccentricity. 
IGoldreich and Tremaind Ill982ll calculated the change in the ec- 
centricity of a ring particle due to a satellite. We are interested 
in the change of the eccentricity of the moonlet induced by a 
ring of particles and therefore swap the satellite mass with the 
mass of a ring particle. Thus, r ewriting their results (Eq. 64 in 
IGoldreich and Tremaind Ill982ll ) in our local notation without 
reference to the semi-major axis, we have 

(25) 



A«5 2 



5.02 m\ G 2 Q 4 b~ 4 . 



Supposing that the moonlet has very small eccentricity, it will 
receive stochastic impulses that cause its eccentricity to undergo 
a random walk that will result in increasing linearly with time, 
so that we may write 



d&\ 



dt 



f 

•J — o> 



W a (b)A& 2 d(l/t b ), (26) 

circulating particles 

where d(l/tb) is the mean encounter rate with particles which 
have an impact parameter in the interval (b, b + db). W a (b) is the 
window function describing the band of particles in circulating 
orbits. Setting d(l/tb) = 3"LCl\b\db/(2m\), we obtain 

dSz 



dt 



= 7.53 | Wa(b)mi\bri<G*Cr 3 db. (27) 

circulating particles 

For high surface densities an additional effect can excite the ec- 
centricity when gravitational wakes occur. The Toomre parame- 
ter Q is defined as 

Q - ^ (28) 

where v is the velocity dispersion of the ring particles^. When the 
surface density is sufficiently high such that Q approaches unity, 



1 The Toomre criterion used here was originally derived for a fiat 
gaseous disc. To make use of it we replace the sound speed by the radial 
velocity dispersion of the ring particles. 
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transient clumps will form in the rings. The typical length scale 
of these structures is given by the critic al Toomre wave length 
A T = 2tt 2 GIQ- 2 (Dais aka et al.L 1200 IB . which can be used to 
estimate a typical mass of the clump: 



= 47r 5 G 2 £ 3 Cr 4 . 



(29) 



Whenever strong clumping occurs, mj should be used in the 
above calculation instead of the mass of a single particle m\\ 



dt 



= 7.53 



circulating clumps 



f 



W' a (b)m T \b\- 3 I.G 2 £l- 3 db, (30) 



where W' a {b) » W a (b) is the appropriate window function. For 
typical parameters used in Sect. [7] this transition occurs at 2 ~ 
200kg/m 2 . 

4.3. Equilibrium eccentricity 

Putting together the results from this and the previous section, 
we can estimate an equilibrium eccentricity of the moonlet. Let 
us assume the eccentricity ez, or the amplitude of the epicyclic 
motion &2, evolves under the influence of excitation and damp- 
ing forces as follows 



dt 



dSj 
dt 

d&l 



, -1 , -1 \p2 

collisions e,horseshoe e,circ / 2 



collisions 



dt 



circulating particles 



(31) 



circulating clumps 



The equilibrium eccentricity is then found by setting the above 
equation equal to zero and solving for &2- 

To make quantitative estimates we need to specify the win- 
dow functions W a (b), Wb+ C (b) and Wd{b) that determine in which 
impact parameter bands the particles are in (see Fig. [I). To com- 
pare our analytic estimates to the numerical simulations pre- 
sented below, we measure the window functions numerically. 
Alternatively, one could simply use sharp cutoffs at some multi- 
ple of the Hill radius. We already made use of this approxima- 
tion as a simple estimate in the previous sections. The results 
may vary slightly, but not significantly. 

However, as the window functions are dimensionless, it is 
possible to obtain the dimensional scaling of £2 by adopting the 
length scale applicable to the impact parameter to be the Hill 
radius r# and simply assume that the window functions are of 
order unity. As already indicated above, all of the circularisation 
times scale as r e oc Q m G l X , or equivalently oc mzKLr^Q). 
Assuming the ring particles have zero velocity dispersion, the 
eccentricity excitation is then due to circulating clumps or parti- 
cles and the scaling of dS^/dt due to this cause is given by Eqs. 
E7]and[30]by 

—± oc mi r H 2 I.G 2 n.- 3 , (32) 

where m, is either m\ or mj. We may then find the scaling of the 
equilibrium value of £2 from consideration of Eqs. [3"Tland[T3"las 



t>2 oc — r 
rri2 



2 

11 ■ 



(33) 



This means that the expected kinetic energy in the non circu- 
lar motion of the moonlet is oc m/r^Q 2 which can be viewed as 



stating that the non circular moonlet motion scales in equiparti- 
tion with the mass m, moving with speed r#£X This speed ap- 
plies when the dispersion velocity associated with these masses 
is zero indicating that the shear across a Hill radius replaces the 
dispersion velocity in that limit. 

In the opposite limit when the dispersion velocity exceeds 
the shear across a Hill radius and the dominant source of eccen- 
tricity excitation is due to collisions, Eq. IBTI gives 



m 2 £2 = mi{&\) 



(34) 



so that the moonlet is in equipartition with the ring particles. 
Results for the two limiting cases can be combined to give an 
expression for the amplitude of the epicyclic motion excited in 
the moonlet of the form 



(35) 



where C, is a constant of order unity. This indicates the tran- 
sition between the shear dominated and the velocity dispersion 
dominated limits. 



5. Processes leading to a random walk in the 
semi-major axis of the moonlet 

We have established estimates for the equilibrium eccentricity 
of the moonlet in the previous section. Here, we estimate the 
random walk of the semi-major axis of the moonlet. In contrast 
to the case of the eccentricity, there is no tendency for the semi- 
major axis to relax to any particular value, so that there are no 
damping terms and the deviation of the semi-major axis from its 
value at time t = grows on average as yft for large t. 

Depending on the surface density, there are different effects 
that dominate the contributions to the random walk of the moon- 
let. For low surface densities, collisions and horseshoe orbits 
are most important. For high surface densities, the random walk 
is dominated by the stochastic gravitational force of circulat- 
ing particles and clumps. In this section, we try to estimate the 
strength of each effect. 

5. 1 . Random walk due to collisions 

Let us assume, without loss of generality, that the guiding centre 
of the epicyclic motion of a ring particle is displaced from the 
orbit of the moonlet by J?li in the inertial frame, whereas in the 
co-rotating frame the moonlet is initially located at the origin 
with = 0. When the impacting particle becomes bound to 
the moonlet, the guiding centre of the epicyclic motion of the 
centre of mass of the combined object, as viewed in the inertial 
frame, is then displaced from the original moonlet orbit by 



- mi m\ 

mi + m2 m.2 



(36) 



This is the analogue of Eq.|22]for the evolution of the semi-major 
axis. For an ensemble of particles with impact parameter b, the 
average centre of epicyclic motion is = b. Thus we can 

write the evolution of due to consecutive encounters as 



dJ{ 2 






dt 


collisions 


\m 2 ) 



mi r°° 3 
ml J-oo 2 



|£| 3 W b+C (b)db, 



(37) 
(38) 



which should be compared to the corresponding expression for 
the eccentricity in Eq.l23l 



6 



H. Rein and J.C.B. Papaloizou: Stochastic orbital migration of moonlets in Saturn's rings 



5.2. Random walk due to stochastic forces from circulating 
particles and clumps 

Particles on a circulating trajectory (see path (a) in Fig. [TJ that 
come close the moonlet will exert a stochastic gravitational 
force. The largest contributions occur for particles within a few 
Hill radii. Thus, we can crudely estimate the magnitude of the 
specific gravitational force (acceleration) due to a single particle 
as 

f = (39) 

fcp {2rH)2 - 

When self-gravity is important, gravitational wakes (or clumps) 
have to be taken into account, as was done in Sec. 14.21 Then a 
rough estimate of the largest specific gravitational forces due to 
self gravitating clumps is 

G TOt 

(2r H + A T ) 

where 2rn has been replaced by 2rn + At- This allows for the 
fact that Aj could be significantly larger than r#, in which case 
the approximate distance of the clump to the moonle t is Aj. 

Following the formalism of iRein and Papaloizou! Il2009ll . we 
define a diffusion coefficient as D — 2r c (/ 2 ), where / is an ac- 
celeration, t c is the correlation time and the angle brackets de- 
note a mean value. We here take the correlation time associated 
with these forces to be the orbital period, 2ji£T 1 . This is the nat- 
ural dynamical timescale of the systems and has been found to 
be a reasonable assumption from an analysis of the simulations 
described below. By considering the rate of change of the en- 
ergy of the moonlet motion, we may estimate the random walk 
in ^ due to circulating particles and sel f-gravitating clumps to 
be given by llRein and Papaloizoul l2Q09tl : 



dJ{ 2 



dt 



dj{ 2 



circulating particles 



4D cp n 



167T Q 



I6n Q 



Gtoi 
(2^F 



- 3 <4> 



dt 



circulating clumps 



\6n Q 



G tot 



(2r H + A T Y 



(41) 



(42) 



(43) 



(44) 



This is only a crude estimate of the random walk undergone 
by the moonlet. In reality several additional effects might also 
play a role. For example, circulating particles and clumps are 
clearly correlated, the gravitational wakes have a large extent 
in the azimuthal direction and particles that spend a long time 
in the vicinity of the moonlet have more complex trajectories. 
Nevertheless, we find that the above estimates are correct up to 
a factor 2 for all the simulations that we performed (see below). 

5.3. Random walk due to particles in horseshoe orbits 

Finally, let us calculate the random walk induced by particles 
on horseshoe orbits. Particles undergoing horseshoe turns on 
opposite sides of the planet produce changes of opposite sign. 
Encounters with the moonlet are stochastic and therefore the 
semi-major axis will undergo a random walk. A single particle 
with impact parameter b will change the semi-major axis of the 
moonlet by 

TO] 



Aj?l = 2- 



TOl + TO 2 



-b. 



(45) 



shear 



main box 



shear 



auxiliary boxes 



Fig. 2. Shearing box, simulating a small patch of a planetary 
ring system. The particles that leave an auxiliary box in the az- 
imuthal (y) direction reenter the same box on the other side and 
get copied into the main box at the corresponding location. All 
auxiliary boxes are equivalent and there are no curvature terms 
present in the shearing sheet. 



Analogous to the analysis in Sec. 15. II the time evolution of J[ is 
then governed by 



dj{ 2 



dt 



horseshoe 



4 S 



W d (b) db. 



(46) 
(47) 



Note that this equation is identical to Eq.[38]except a factor 4, as 
particles with impact parameter b will leave the vicinity of the 
moonlet at -b. 



6. Numerical Simulations 

We perform realistic three dimensional simulations of Saturn's 
rings with an embedded moonlet and verify the analytic esti- 
mates presented above. The nomenclature and parameters used 
for the simulations are listed in Table Q] 



6.1. Methods 

The gravitational force s are calculated with a Barnes-Hut tree 
llBarnes and Hutt 1986ll. T he numerical scheme is similar to that 
used by Rein et al. Il2010ll . Additionall y, we implemented a sym- 
plectic integrator for Hill's equations flOuinn et all I2010I1 . This 
turned out to be beneficial for accurate energy conservation at 
almost no additional cost when the eccentricity of the moonlet 
(~ 1CT 8 ) was small and integrations were undertaken over many 
hundreds of dynamical timescales. 

To further speed up the calculations, we run two coupled 
simulations in parallel. The first adopts the main box which in- 
corporates a moonlet. The second adopts an auxiliary box which 
initially has the same number of particles but which does not 
contain a moonlet. This is taken to be representative of the unper- 
turbed background ring. We describe this setup as enabling us to 
adopt pseudo shear periodic boundary conditions. We consider 
the main box together with eight equivalent auxiliary boxes to 
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Name 


S 




dt 






X Ly 




N 


EQ5825 


50 kg/m 2 


Z3 


m 


4s 


1000 


m 


X 


1 ADA 


m 


7.2k 


EQ5858 


50kg/m 2 


^ a 


m 


4s 


1000 


m 


X 


1 AAA 
1UU0 


m 


7.2k 


EQ28825 


200 kg/m 2 


£J 


m 


4s 


1000 


m 


X 


1 AAA 
1UU0 


m 


28.8k 


EQ28858 


200 kg/m 2 


50 


m 


4s 


1000 


m 


X 


1000 


m 


28.8k 


EQ48825 


400 kg/m 2 


25 


m 


4s 


1000 


m 


X 


1000 


m 


57.6k 


EQ48858 


400 kg/m 2 


50 


m 


4s 


1000 


m 


X 


1000 


m 


57.6k 


EQ48858DT 


400 kg/m 2 


50 


m 


40s 


1000 


m 


X 


1000 


m 


57.6k 


EQ5858DTW 


50 kg/m 2 


50 


m 


40s 


2000 


m 


X 


2000 


m 


28.8k 


EQ28858DTW 


200 kg/m 2 


50 


m 


40s 


2000 


m 


X 


2000 


m 


115.2k 


EQ48858DTW 


400 kg/m 2 


50 


m 


40s 


2000 


m 


X 


2000 


in 


230.0k 



Table 1. Initial simulation parameters. The second column gives the surface density of the ring. The third gives the moonlet radius. 
The fourth column gives the time step. The fifth and sixth columns give the lengths of the main box as measured in the .ry-plane and 
the number of particles used, respectively. 



be stacked as illustrated in Fig. [2] All auxiliary boxes are identi- 
cal copies whose centres are shifted according to the background 
shear as in a normal shearing sheet. On account of this motion 
auxiliary boxes are removed when their centres are shifted in 
azimuth by more than l.5L y from the centre of the main box 
and then reinserted in the same orbit on the opposite side of the 
main box so that the domain under consideration is prevented 
from shearing out. If a particle in the main box crosses one of 
its boundaries, it is discarded. If a particle in the auxiliary box 
crosses one of its boundaries, it is reinserted on the other side 
of this box, according to normal shear periodic boundary con- 
ditions. But in addition it is also copied into the corresponding 
location in the main box. We describe this procedure as applying 
pseudo shear periodic boundary conditions. 

In a similar calculation. Lewis and Steward Il2009ll use a very 
long box (about 10 times longer than the boxes used in this pa- 
per) to ensure that particles are completely randomised between 
encounters with the moonlet. We are not interested in the long 
wavelength response that is created by the moonlet. Effects that 
are most important for the moonlet's dynamical evolution are 
found to happen within a few Hill radii. Using the pseudo shear 
periodic boundary conditions, we ensure that incoming particles 
are uncorrected and do not contain prior information about the 
perturber. This setup speeds up our calculations by more than an 
order of magnitude. 

The gravity acting on a particle in the main box, which also 
contains the moonlet, is calculated by summing over the parti- 
cles in the main box and all auxiliary boxes. The gravity acting 
on a particle in the auxiliary box is calculated the standard way, 
by using ghost boxes which are identical copies of the auxiliary 
box. Tests have indicated that our procedure does not introduce 
unwanted fluctuations in the gravitational forces. 

The moonlet is allowed to move freely in the main box. 
However, in order to prevent it leaving the computational do- 
main, as soon as the moonlet has left the innermost part (de- 
fined as extending one eighth of the box size), all particles are 
shifted together with the box boundaries, such that the moonlet 
is returned to the centre of the box. This is possible because the 
shearing box approximation is invariant with respect to transla- 
tions in the xy plane (see Eq.[TJi. 

Collisions between particles are resolved using the instan- 
taneous collision m odel and a velocity d ependent coefficient of 
restitution given by [Bridges et al.llll984ll : 




e(v) = min 



0.34 ■ 



lcm/s 



-0.234 



. 1 



(48) 



Fig. 3. Snapshots of the particle distribution and the moonlet 
(blue) in simulation EQ400 5 0DTW. The wake is much clearer than 
in Fig. |4]as the box size is twice as large. 



where v\\ is the impact speed projected on the axis between the 
two particles. The already ex isting tree struct ure is reused to 
search for nearest neighbours llRein et all 1201 Oil . 



6.2. Initial conditions and tests 

Throughout this paper, we use a distribution of particle sizes, r\ , 
ranging from lm to 5m with a slope of q — -3. Thus dN/dri oc 
r v q . The density of both the ring particles and the moonlet is 
taken to be 0.4 g /cm 2 . This is at the lower e nd of what has been 
assum ed reasonable for Saturn's A ring llLewis and StewartL 
l2009ll . The moonlet radius is taken to be either 50 m or 25 m. 
We found that using a larger ring particle and moonlet density 
only leads to more particles being bound to the moonlet. This 
effectively increases the mass of the moonlet (or equivalently its 
Hill radius) and can therefore be easily scaled to the formalism 
presented here. Simulating a gravitational aggregate of this kind 
is computationally very expensive, as many more collisions have 
to be resolved each time-step. 
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(c) S = 400 kg/m 2 , r 2 = 25 m (d) S = 400 kg/m 2 , r 2 = 50 m 

Fig. 4. Snapshots of the particle distribution and the moonlet (blue) for different surface densities after 25 orbits. The wake created 
by the moonlet is much longer than the size of the box and therefore hardly visible in these images. 
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The initial velocity dispersion of the particles is set to 
lmm/s for the x and y components and 0.4mm/s for the z 
component. The moonlet is placed at a semi-major axis of 
a = 130000 km, corresponding to an orbital period of P — 
13.3 hours. Initially, the moonlet is placed in the centre of the 
shearing box on a circular orbit. 

We run simulations with a variety of moonlet sizes, sur- 
face densities and box sizes. The dimensions of each box, as 
viewed in the (x, y) plane, are specified in Table Q] Because of 
the small dispersion velocities, the vertical motion is automati- 
cally strongly confined to the mid-plane. After a few orbits, the 
simulations reach an equilibrium state in which the velocity dis- 
persion of particles does not change anymore. 

We ran several tests to ensure that our results are converged. 
Simulation EQ40858DT uses a ten times larger time step than 
simulation EQ4QQ5Q. Simulations EQ5Q58DTW, EQ2QQ5QDTWand 
EQ4QS5SDTWuse a box that is twice the size of that used in simu- 
lation EQ4QQ5Q. Therefore, four times more particles have been 
used. No differences in the equilibrium state of the ring particles, 
nor in the equilibrium state of the moonlet have been observed 
in any of those test cases. 

A snapshot of the particle distribution in simulation 
EQ4QQ5QDTW is shown in Fig. [3] Snapshots of simulations 
EQ20Q25, EQ2QQ5Q, EQ4QQ25 and EQ4QQ5Q, which use the 
smaller box size, are shown in Fig. [4] In all these cases, the 
length of the wake that is created by the moonlet is much 
longer than the box size. Nevertheless, the pseudo shear peri- 
odic boundary conditions allow us to simulate the evolution of 
the moonlet accurately. 

7. Results 

7.1. Impact bands 

The impact band window functions W„{b), Wb+ C (b) and Wd(b) 
are necessary to estimate the damping timescales and the 
strength of the excitation. To measure those, each particle that 
enters the box is labelled with it's impact parameter. The possi- 
ble outcomes are plotted in Fig.Q] (a) being a circulating parti- 
cle which leaves the box on the opposite side, (b) representing a 
collision with other ring particles close to the moonlet where the 
maximum distance from the centre of the moonlet has been taken 
to be twice the moonlet radius, (c) being a direct collision with 
the moonlet and finally (d) showing a horseshoe orbit in which 
the particle leaves the box on the same side that it entered. 

The impact band (b) is considered in addition to the impact 
band (c) because some ring particles will collide with other ring 
particles that are (temporarily) bound to the moonlet. Thus, these 
collisions take part in the transfer the energy and momentum 
to the moonlet. Actually, if the moonlet is simply a rubble pile 
of ring particles as suggested by iPorco et al.l ll2007ll . then there 
might be no solid moonlet core and all collisions are in the im- 
pact band (b). 

Fig. [5] shows the impact bands, normalised to the Hill ra- 
dius of the moonlet rn for simulations EQ2S050 and EQ200250. 
For comparison, we also plot sharp cutoffs, at 1.5r# and 2.5r#. 
Our results show no sharp discontinuity because of the velocity 
dispersion in the ring particles which is ~ 5 mm/s. This corre- 
sponds to an epicyclic motion of ~ 40 m or equivalently ~ 0.8 r# 
and ~ 1 .6 r# (for rn = 50m and r> = 25m, respectively), which 

2 Y,i=abcd Wj > 1 because all particles are shifted from time to time 
when the moonlet is to far away from the origin (see Sect. l6.lt . In this 
process, particles with the same initial impact parameter might become 
associated with the more than one impact band. 
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Name 


Numerical results 


Analytic results 






collisions 


horseshoe 


EQ5025 


18.9 


18.4 


1712 


EQ5058 


27.6 


33.8 


3424 


EQ20825 


9.0 


3.9 


428 


EQ2085S 


12.9 


8.5 


856 


EQ48025 


2.8 


2.1 


214 


EQ40050 


5.5 


4.4 


428 



Table 2. Eccentricity damping timescale r e of the moonlet in 
units of the orbital period. The second column lists the simula- 
tion results. The third and fourth column list the analytic esti- 
mates of collisional and horseshoe damping timescales, respec- 
tively. 



explains the cutoff width of approximately one Hill radii. The 
impact bands are almost independent of the surface density and 
depend only on the Hill radii and the mean epicyclic amplitude 
of the ring particles, or, more precisely, the ratio thereof. 

7.2. Eccentricity damping timescale 

To measure the eccentricity damping timescale, we first let the 
ring particles and the moonlet reach an equilibrium and integrate 
them for 200 orbits. We then change the velocity of the moonlet 
and the ring particles within 2ru- The new velocity corresponds 
to an eccentricity of 6 • 10~ 7 , which is well above the equilib- 
rium value. We then measure the decay timescale r e by fitting a 
function of the form 

e d (t) = (e eq ) + (6 ■ 10- 7 - <«„>) e- ,lT '. (49) 

The results are given in Table|2] They are in good agreement with 
the estimated damping timescales from Sect. 13. Il and l3.2[ show- 
ing clearly that the most important damping process, in every 
simulation considered here, is indeed through collisional damp- 
ing as predicted by comparing Eq. [15] with Eqs.[T9land[2T1 

7.3. Mean moonlet eccentricity 

The mean eccentricity of the moonlet is measured in all simula- 
tions after several orbits when the ring particles and the moonlet 
have reached an equilibrium state. To compare this value with 
the estimates from Sect. [2] we set Eq. [3T1 equal to zero and use 
the analytic damping timescale listed in Table|2] The analytic es- 
timates of the equilibrium eccentricity are calculated for each ex- 
citation mechanism separately to disentangle their effects. They 
are listed in the third, fourth and fifth column in Table [3] The 
sixth column lists the analytic estimate for the mean eccentricity 
using the sum of all excitation mechanisms. 

For all simulations, the estimates are correct within a factor 
of about 2. For low surface densities, the excitation is dominated 
by individual particle collisions. For larger surface densities, it 
is dominated by the excitation due to circulating self gravitat- 
ing clumps (gravitational wakes). The estimates and their trends 
are surprisingly accurate, as we have ignored several effects (see 
below). 

7.4. Random walk in semi-major axis 

The random walk of the semi-major axis a (or equivalently the 
centre of epicyclic motion J[ in the shearing sheet) of the moon- 
let in the numerical simulations are measured and compared to 
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r 2 =50m r 2 =25m 



o 




impact parameter [r HiN ] impact parameter [r HN |] 



Fig. 5. Impact bands for simulation EQ2005Q (left) and EQ20025 (right). The vertical lines indicate sharp cutoffs at 1.5r# and 2.5rn- 



the analytic estimates presented in Sect. [3] We ran one simula- 
tion per parameter set. To get a statistically meaningful expres- 
sion for the average random walk after a given time yi(Af), we 
average all pairs of J[{t) and ^l(f') for which t-t' - At. In other 
words, we assume the system satisfies the Ergodic hypothesis. 

J{(At) then grows like VAf and we can fit a simple square 
root function. This allows us to accurately measure the average 
growth in J[ after Af = 100 orbits by running one simulation for 
a long time and averaging over time, rather than running many 
simulations and performing an ensemble average. The measured 
values are given in the second column of Table [4] 

The values that correspond to the analytic expressions in 
Eqs. [47] [38] @2] and [44] are listed in columns three, four, five 
and six, respectively. 

For low surface densities, the evolution of the random walk 
is dominated by collisions and the effect of particles on horse- 
shoe orbits. For large surface densities, the main effect comes 
from the stochastic gravitational force due to circulating clumps 
(gravitational wakes). 

7.5. Long term evolution and observability 

The change in semi-major axis is relatively small, a few tens of 
meters after 100 orbits (= 50 days). The change in longitude 
that corresponds to this is, h owever, is much larger. We there- 
fore extend the discussion of iRein and Papaloizoul Il2009ll . who 
derive the time evolution of orbital parameters undergoing a ran- 
dom walk, to the time evolution of the mean longitude. This is 
of special interest in the current situation, as in observations of 
Saturn's rings, it is much easier to measure the shift in the longi- 
tude of an embedded moonlet than the change in its semi-major 
axis. 

The time derivative of the mean motion is given by 



3 

n — — 

2 \ a 



—a - - — 



a 



(50) 



200 



150 



100 



r so 



-50 



-100 



-150 



-200 



SF*50kg/m* r 2 =25m 
.'■£=50kg/rrC r 2 =50m 
/L=200kg/rrC r 2 =25m 
/ L=200kg/rrC r 2 =50m 
L=400kg/rrC r 2 =25m 
.- - -E=400kg/m r 2 =50m 



0.1 



0.2 0.3 
time [years] 



0.4 



0.5 



Fig. 6. Offset in the azimuthal distance due to the random walk in 
the semi-major axis a measured in simulations EQ5025, EQ505Q, 
EQ2Q025, EQ2QQ58, EQ4Q025 and EQ40050. Individual moon- 
lets may show linear, constant or oscillatory growth. On average, 
the azimuthal offset grows like f 3 ^ 2 . 



the double integral 

A(t) = f (no + An(t')) dt' 
Jo 



not _ rr"w w 

Jo Jo a 



(51) 
(52) 



The root mean square of the difference compared to the unper- 
turbed orbit (no) is 



(AA) 2 = ((A(t)-nt) 2 ) 

'' ' 9F g (t")F e (t"") 



- ffff 



dt"" dt'" dt" dt' 



(53) 
(54) 



where Fg is the time dependent stochastic force in the 8 (or in 
shearing sheet coordinates, y) direction and we have assumed a 
nearly circular orbit, e«l, The mean longitude is thus given by 



t,t',t,t" 



g(\t" -t""\)dt""dt"'dt"dt' (55) 
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Name 


Numerical results 


collisions 


Analytic results 
circulating particles circulating clumps 


total 


EQ5025 


6.1 ■ 1(T 8 


4.5 • 1(T 8 


1.9 • 1(T 8 


0.3 • 10- 8 


4.9 • 10 


■8 


EQ5050 


4.3 • 1(T 8 


1.6 • lfr 8 


1.4 ■ 10 8 


0.2 • 10- 8 


2.1 • 10 


-8 


EQ20025 


6.4 • irr 8 


4.6 • 10~ 8 


1.7 ■ lfr 8 


2.0 • 10~ 8 


5.3 • 10 


-8 


EQ20050 


4.9 • 10~ 8 


1.6 - 10 s 


1.4 ■ 10 8 


1.6 • 10" 8 


2.7 • 10 


-8 


EQ40025 


n.3 • irr 8 


4.6 • lfr 8 


2.5 • lfr 8 


8.4 • 10- 8 


9.9 • 10 


-8 


EQ40050 


7.2 • 10~ 8 


1.6- 1(T 8 


1.5 ■ lfr 8 


4.9 • 10~ 8 


5.4- 10 


-8 



Table 3. Equilibrium eccentricity (e e q) of the moonlet. The second column lists the simulation results. The third, fourth and fifth 
column list the analytic estimates of the equilibrium eccentricity assuming a single excitation mechanism. The last column lists the 
analytic estimate of the equilibrium eccentricity summing over all excitation mechanisms. 



Name 


Numerical results 


horseshoe 


collisions 


Analytic results 
circulating particles 


circulating clumps 


total 


EQ5025 


22.7m 


9.0m 


10.5m 


17.6m 


0.3m 


22.4m 


EQ5050 


12.6m 


4.5m 


4.8m 


4.4m 


0.1m 


7.9m 


EQ20025 


38.1m 


18.1m 


25.1m 


17.6m 


15.7m 


38.9m 


EQ20050 


22.7m 


13.6m 


9.6m 


4.4m 


4.8m 


14.7m 


EQ40025 


78.1m 


25.6m 


36.5m 


17.6m 


86.6m 


100.7m 


EQ40050 


45.6m 


12.8m 


13.4m 


4.4m 


31.4m 


36.7m 



Table 4. Change in semi-major axis of the moonlet after 100 orbits, J[(At = 100 orbits). The second column lists the simula- 
tion results. The third till sixth columns list the analytic estimates of the change in semi-major axis assuming a single excitation 
mechanism. The last column lists the total estimated change in semi-major axis summing over all excitation mechanisms. 



Note that Eq.[62]is an exact solution, whereas Eq.[57]is a statisti- 
cal quantity, describing the root mean square value in an ensem- 
ble average. In the stochastic and laminar case, (A/1) 2 grows like 
~ f 3 and ~ f 4 , respectively. This behaviour allows to discrim- 
inate between a stochastic and laminar migration by observing 
AA over an extended period of time. 

Results from individual simulations (i.e. not an ensemble av- 
erage) are plotted in Fig. [6] Here, AA is expressed in terms of the 
azimuthal offset relative to a Keplerian orbit. One can see that 
for an individual moonlet, the shift in azimuth can appear to be 
linear, constant or oscillating on short timescales (see curve for 
simulation EQ2Q025). This is partly because of the lower order 
terms in Eq. [56] Though, on average, the azimuthal offset grows 
very rapidly, as ~ f 3 ' 2 (see Eq. [57l i for t » r. 

8. Conclusions 

In this paper, we have discussed the dynamical response of an 
embedded moonlet in Saturn's rings to interactions with ring 
particles both analytically and by the use of realistic three dimen- 
sional many particle simulations. Both the moonlet and the ring 
particle density were taken to be 0.4g/cm 2 . Moonlets of radius 
25 m and 50 m were considered. Particle sizes ranging between 
1 m and 5 m were adopted. 

We estimated the eccentricity damping timescale of the 
moonlet due to collisions with ring particles and due to the 
response of ring particles to gravitational perturbations by the 
moonlet analytically. We found the effects due to the response of 
particles on horseshoe and circulating orbits are negligible. On 
the other hand the stochastic impulses applied to the moonlet by 
circulating particles were found to cause the square of the eccen- 
tricity to grow linearly with time as did collisions with particles 
with non zero velocity dispersion. A balance between excitation 
and damping processes then leads to an equilibrium moonlet ec- 
centricity. 

We also estimated the magnitude of the random walk in the 
semi-major axis of the moonlet induced by collisions with indi- 
vidual ring particles and the gravitational interaction with parti- 



= 9 J&L (-2t 4 + (2r 3 f + 2r 4 + r¥) + irf 3 ) (56) 

where g(t) is the auto-correlation function of Fg, which has been 
assumed to be exponentially decaying with an auto-correlation 
time t, thus g(f) = exp(-[f|/r). In the limit t » r, the above 
equation simplifies to 



3«> 



T . 



(AA) 2 = t 3 = (An) 2 ■-. (57) 

a 1 6 

On shorter timescale, the other terms in Eq. [56] are important, 
leading to linear and oscillatory behaviour. For uncorrected 
noise (e.g. collisions) Eq.[57]is replaced by 



(A/1) 2 = ^ P/r, (58) 



where (Av) is the average velocity change per impulsive event 
and t is taken to be the average time between consecutive events. 

Instead of a stochastic migration, let us also assume a lami- 
nar migration with constant migration rate r a , defined by 

a 3 n 

T„ = 7-7- (59) 

a In 

As above, we can calculate the longitude as a function of time as 
the following double integral 



Jo 



/1(f) = I (n + An(f')) dt' (60) 

n'' 3 n 3 n 

- — dt"dt' = n t - -— t 2 (61) 
2 T a 4 T a 

and thus 

(A/1) 2 = t 4 = (A«) 2 ^. (62) 

16 Tt 4 
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cles and gravitational wakes. There is no tendency for the semi- 
major axis to relax to any particular value, so that there are no 
damping terms. The deviation of the semi-major axis from its 
value at time t — grows on average as yfi for large t . 

From our simulations we find that the evolution of the eccen- 
tricity is indeed dominated by collisions with ring particles. For 
large surface densities (more than 200kg/m 2 ) the effect of grav- 
itational wakes also becomes important, leading to an increase 
in the mean steady state eccentricity of the moonlet. When the 
particle velocity dispersion is large compared to Q.rn, we obtain 
approximate energy equipartition between the moonlet and ring 
particles as far as epicyclic motion is concerned. 

Similarly, the random walk in the semi-major axis was found 
to be dominated by collisions for low surface densities and by 
gravitational wakes for large surface densities. In addition we 
have shown that on average, the difference in longitude AA of 
a stochastically forced moonlet grows with time like P^ 2 for 
large t. 

The distance travelled within 100 orbits (50 days at a dis- 
tance of 130000km) is, depending on the precise parameters, 
of the order of 10- 100m. This translates to a shift in longitude 
of several hundred kilometres. The shift in A is not necessarily 
monotonic on short timescales (see Fig. [6]). We expect that such a 
shift sh ould be easily detectab le by the Cassini spacecraft. And 
indeed, iTiscareno et al.l Il2010ll report such an observation, al- 
though the data is not publicly available yet, at the time when 
this paper was submitted. 

There are several effects that have not been included in this 
study. The analytic discussions in Sec. [3] Sec.|4]and Sec.[5]do not 
model the motion of particles correctly when their trajectories 
take them very close to the moonlet. Material that is temporarily 
or permanently bound to the moonlet is ignored. The density 
of both the ring particles and the m oonlet are at the lower en d 
of what is assumed to be reasonable IlLewis and Steward l2009ll . 
This has the advantage that for our computational setup only a 
few particles (with a total mass of less than 10% of the moonlet) 
are bound to the moonlet at any given time. In a future study, 
we will extend the discussion presented in this paper to a larger 
variety of particle sizes, moonlet sizes and densities. 
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Appendix A: Response calculation of particles on 
horseshoe orbits 

A. 1 . Interaction potential 

Due to some finite eccentricity, the moonlet undergoes a small 
oscillation about the origin. Its Cartesian coordinates then be- 
come (X, Y), with these being considered small in magnitude. 
The components of the equation of motion for a ring particle 
with Cartesian coordinates (x,y) = ri are 

dt 1 dt m\ ox 



and 



(A.l) 



^ +2Q ^ = _J_^1£, (A . 2) 

dt 1 dt m\ oy 

where the interaction gravitational potential due to the moonlet 
is 



1.2 



Gm\rri2 



with the cylindrical coordinates of the particle and moonlet be- 
ing (r, 0) and (R, O) respectively. This may be expanded correct 
to first order in R/r in the form 



T-1,2 = - 



Gm\ni2 Gm\ miR cos(0 — <D) 



(A.4) 



The moonlet undergoes small amplitude epicyclic oscillations 
such that X - &2 cos(£2? + e), Y = -2&2 sin(Q? + e) where e is its 
small eccentricity and e is an arbitrary phase. Then W\ t 2 may be 
written as 



1.2 



where 



vp' 

M,2 



1,2' 



Gm\mir &i(x cos(£2? + e) - 2y sin(Qf + e)) 



(A.5) 



(A.6) 



gives the part of the lowest order interaction potential associated 
with the eccentricity of the moonlet. 

We here view the interaction of a ring particle with the moon- 
let as involving two components. The first, due to the first term 
on the right hand side of Eq. IA. 41 operates when £2=0 and re- 
sults in standard horseshoe orbits for ring particles induced by 
a moonlet in circular orbit. The second term, Eq. IA.6I perturbs 
this motion when £2 is small. We now consider the response of 
a ring particle undergoing horseshoe motion to this perturbation. 
In doing so we make the approximation that the variation of the 
leading order potential Eq. IA.4l due to the induced ring particle 
perturbations may be neglected. This is justified by the fact that 
the response induces epicyclic oscillations of the particle which 
are governed by the dominant central potential. 

A.2. Response calculation 

Setting x — > x + £ t , y — > y + w here f is the small response 
displacement induced by Eq. I A. 61 and linearising Eqs. IA.2I we 
obtain the following equations for the components of £ 



dt 1 dt m\ ox 



^| + 2Q^ 
dt 2 dt 



1 W'v 



m\ dy 

From these we find a single equation for % x in the form 
-fit 1 trv.. r 1 &V[ 2 



df- 



ni\ d. 



m\ dy 



dt = F. 



{A.l) 
(A.8) 

(A.9) 



When performing the time integral on the right hand side of the 
above, as we are concerned with a potentially resonant epicyclic 
response, we retain only the oscillating part. 

The solution to Eq. lA.9l which is such that £ vanishes in the 
distant past (f = -00) when the particle is far from the moonlet 
may be written as 



£t = acos(Qf) +/?sin(Qf), 
where 

Fsin(Qf) 



(r 2 + R 2 - 2rRcos(cf> - (D)) 1 / 2 



(A.3) p 



-f 

U -c 



Q 

Fcos(Qf) 



dt and 



dt. 



(A. 10) 

(A.ll) 
(A. 12) 
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After the particle has had its closest approach to the moonlet 
and moves to a large distance from it it will have an epicyclic 
oscillation with amplitude and phase determined by 



U —i 



Ax 



f 



x F sin(Of) 
Fcos(Qf) 



dt and 



Q 



dt. 



(A.13) 
(A. 14) 



In evaluating the above we note that, although F vanishes when 
the particle is distant from the moonlet at large \t\, it also oscil- 
lates with the epicyclic angular frequency O which results in a 
definite non zero contribution. This is the action of the co-orbital 
resonance. It is amplified by the fact that the encounter of the 
particle with the moonlet in general occurs on a horseshoe libra- 
tion time scale which is much longer than Q _I . We now consider 
this unperturbed motion of the ring particles. 

A3. Unperturbed horseshoe motion 

The equations governing the unperturbed horseshoe motion are 
Bqs.rOwith 



^1.2 = <2 = " 



(A.15) 



We assume that this motion is such that x varies on a time scale 
much longer than Q. 1 so that we may approximate the first of 
these equations as x = -2/(3Q.)(dy/df). Consistent with this we 
also neglect x in comparison to y in Eq. IA. 151 From the second 
equation we than find that 



d*y = 3 fflf?, 2 
dt 2 mi dy 

which has a first integral that may be written 

tdW _ 6Gm 2 9Q. 2 b 2 
\dt) " \yT + 4 ' 



(A. 16) 



(A. 17) 



where as before b is the impact parameter, or the constant 
value of \x\ at large distances from the moonlet. The value of 
y for which the horseshoe turns is then given by y — yo = 
24Gm2/(9Qrb 2 ). At this point we comment that we are free to 
choose the origin of time t — such that it coincides with the 
closest approach at which y - yo- Then in the approximation we 
have adopted, the horseshoe motion is such that y is a symmetric 
function of t. 

A.4. Evaluation of the induced epicyclic amplitude 



Using Eqs. IA.6I IA.9I and IA. 171 we may evaluate the epicyclic 
amplitudes from Eq. IA. 14-1 In particular we find 



Gm 2 6 2 cos(Qf + e) / (3x 2 + I2y 2 ) 



(A. 18) 



When evaluating Eq. IA.141 consistent with our assumption that 
the epicyclic oscillations are fast, we average over an orbital pe- 
riod assuming that other quantities in the integrands are fixed, 
and use Eq. IA. 171 to express the integrals with respect to t as 
integrals with respect to y. We then find = -A <x, sin e and 
y6oo = -A a, cos e where 



A 



-r 



^6Gni2yo &i 

6J> 3 VI -yo/y 



dy, 



(A. 19) 



with r 2 = 8Gm2/(3Q 2 yo)(l —yo/y) +y 2 - From this we may write 
Gmi So 

|A»l = Si/ = ^rm, (A.20) 



where &\ / is the final epicyclic motion of the ring particle and 
the dimensionless integral I is given by 



I = 



Jyn 



2 -3 

y r 



6Gm 2 J yo yi -y Q /y 



dy. 



(A.21) 



We remark that the dimensionless quantity rj = %Gni2l{'i^l 2 y i Q ) = 
&( r H/yo) 3 ', with ru being the Hill radius of the moonlet. It is re- 
lated to the impact parameter b by rj — 2~ 6 (b/rn) 6 '. Thus for an 
impact parameter amounting to a few Hill radii, in a very ap- 
proximate sense, 77 is of order unity, yo is of order r# and the 
induced eccentricity £1 y is of order &2- 
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